In this script I’ll be looking at how to analyse the outputs from the MCMCglmm models to get the something about the elaboration/exploration stories.
I’m going to look at it using two approaches: * the blob-wise one where I’ll look at the differences between whole blobs (ellipses) in terms of elaboration/exploration (Robinson & Beckerman style); * and the tip-wise one where I’ll look at the elaboration/exploration score of the tips on the tree (cf. their blobs) relative to the whole phylogeny and to their respective blobs (clades).
I’m going to focus on the data from Gavin and especially the models 5 and 6 that are:
model_phylo1_clade3 with one phylogenetic random effect (across the whole tree) and three clade residual effects.model_phylo4 with three phylogenetic random effect (one for each clade) and three clade residual effects.## Loading the correct models
load("../Data/Processed/model_list.rda")
model_phylo1_clade3 <- model_list[[4]]
model_phylo4 <- model_list[[7]]
And here’s what the trait space looks like:
## Loading the data
load("../Data/Processed/morphdat.rda")
## Plotting it
colour_vector <- c("orange", "blue", "darkgreen")
plot(morphdat[, c(1,2)], pch = 19,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
legend("topleft", legend = levels(morphdat$clade), col = colour_vector, pch = 19)
Note that the PC% are from the whole PC in Gavin’s example.
Here we’re basically comparing multidimensional ellipses (here 3D ones but I think we should push it to more dimensions!).
First we can visualise some of these ellipses (100, randomly drawn from the posterior):
source("../Functions/plot.ellipses.R")
source("../Functions/get.covar.R")
## Selecting 100 random covariance matrices from the MCMCglmm
covar_matrices <- get.covar(model_phylo1_clade3,
n = 100,
simplify = FALSE)
## Plotting the results
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "none",
col = c("grey",colour_vector),
add = TRUE)
Because this is really messy, we can centre these ellipses on the groups ellipses average centres:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "level",
col = c("grey",colour_vector),
add = TRUE)
We can then measure different stuff from these ellipses (e.g. going full statistics from Robinson & Beckerman).
However, here I’m just going to focus on the elaboration/exploration of these ellipses. For the three clades, I’m going to measure:
For these three metrics I expect:
Note that the three metrics are measured independently of the ellipse position in space and are measured in 3D.
source("../Functions/get.axes.R")
## Get all the phylo main axes
level_centred_major_axes <- get.axes(covar_matrices, centre = "level")
uncentered_major_axes <- get.axes(covar_matrices, centre = "none")
And this is what it looks like for the uncentred ones:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "none",
col = c("grey",colour_vector),
add = TRUE)
plot.all.axes(uncentered_major_axes,
col = c("grey",colour_vector),
add = TRUE)
And for the centred ones:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(covar_matrices, centre = "level",
col = c("grey",colour_vector),
add = TRUE)
plot.all.axes(level_centred_major_axes,
col = c("grey",colour_vector),
add = TRUE)
Note that the major axes look slightly off for some of the ellipses. I think this is due to the difference in aspect ratio in both plots that exagerates vertical positions and deforms the ellipses:
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)",
main = "Major axes of the ellipses\n(aspect ratio = 1)",xlim = c(-1.5, 1.5), ylim = c(-1.5, 1.5))
plot.ellipses(covar_matrices, centre = "level",
col = c("grey",colour_vector),
add = TRUE)
plot.all.axes(level_centred_major_axes,
col = c("grey",colour_vector),
add = TRUE)
This is an example with just the major axes (projected in 2D):
plot.all.axes(uncentered_major_axes,
col = c("grey",colour_vector),
add = FALSE, main = "Major axes (un-centred)")
We can then compare these vectors per blobs. We in order to facilitate the vectors comparisons, we will translate the groups vectors (clade ones) on the origin of the phylo vector.
source("../Functions/vector.diff.R")
## Get all the angles and stuff
gull_results <- mapply(vector.diff, uncentered_major_axes$cladegulls, uncentered_major_axes$phylo)
plovers_results <- mapply(vector.diff, uncentered_major_axes$cladeplovers, uncentered_major_axes$phylo)
sandpipers_results <- mapply(vector.diff, uncentered_major_axes$cladesandpipers, uncentered_major_axes$phylo)
## Sort the results per type
angles <- cbind(gull_results["angle", ],
plovers_results["angle", ],
sandpipers_results["angle", ])
projections <- cbind(gull_results["projection", ],
plovers_results["projection", ],
sandpipers_results["projection", ])
rejections <- cbind(gull_results["rejection", ],
plovers_results["rejection", ],
sandpipers_results["rejection", ])
colnames(angles) <- colnames(projections) <- colnames(rejections) <- levels(morphdat$clade)
## Plot the blob wise metrics
par(mfrow = c(3, 1))
boxplot(angles, col = colour_vector, ylab = "angle")
boxplot(projections, col = colour_vector, ylab = "projection (scaled)")
boxplot(rejections, col = colour_vector, ylab = "rejection (scaled)")
Same but with absolute values and angles modulo 90:
## Plot the blob wise metrics
par(mfrow = c(3, 1))
boxplot(angles %% 90, col = colour_vector, ylab = "angle (modulo 90)")
boxplot(abs(projections), col = colour_vector, ylab = "projection (scaled - absolute)")
boxplot(abs(rejections), col = colour_vector, ylab = "rejection (scaled - absolute)")
For this example, the MCMCglmm model has four random (phylo) effects (one for the whole phylo and one per clade) rather than just the one for the whole phylogeny. The model seems to struggle converging and there are things off with the centring (they seem to be all centred no matter what) but that last bit should not be a problem for comparing the vectors since they are translated anyways.
plot(morphdat[, c(1,2)], pch = 19, cex = 0.5,
col = colour_vector[morphdat$clade],
xlab = "PC1 (90.5%)",
ylab = "PC1 (6.86%)")
plot.ellipses(model_phylo4, n = 100,
centre = "none", add = TRUE, col = c("grey",colour_vector, "pink"))
The pink ellipse is representing the residuals.
And we can still do the whole axis stuff, same way as above:
## Get the covariance matrices
covar_matrices2 <- get.covar(model_phylo4,
n = 100,
simplify = FALSE)
## Drop the residuals
covar_matrices2 <- covar_matrices2[-5, ]
## Get all the phylo main axes
uncentered_major_axes2 <- get.axes(covar_matrices2, centre = "none")
names(uncentered_major_axes2) <- c("phylo", "cladegulls", "cladeplovers", "cladesandpipers")
Here’s what the axes look like without the ellipses and stuff:
plot.all.axes(uncentered_major_axes2,
col = c("grey",colour_vector),
add = FALSE, main = "Major axes (centred)")
source("../Functions/vector.diff.R")
## Get all the angles and stuff
gull_results <- mapply(vector.diff, uncentered_major_axes2$cladegulls, uncentered_major_axes2$phylo)
plovers_results <- mapply(vector.diff, uncentered_major_axes2$cladeplovers, uncentered_major_axes2$phylo)
sandpipers_results <- mapply(vector.diff, uncentered_major_axes2$cladesandpipers, uncentered_major_axes2$phylo)
## Sort the results per type
angles <- cbind(gull_results["angle", ],
plovers_results["angle", ],
sandpipers_results["angle", ])
projections <- cbind(gull_results["projection", ],
plovers_results["projection", ],
sandpipers_results["projection", ])
rejections <- cbind(gull_results["rejection", ],
plovers_results["rejection", ],
sandpipers_results["rejection", ])
colnames(angles) <- colnames(projections) <- colnames(rejections) <- levels(morphdat$clade)
## Plot the blob wise metrics
par(mfrow = c(3, 1))
boxplot(angles, col = colour_vector, ylab = "angle")
boxplot(projections, col = colour_vector, ylab = "projection (scaled)")
boxplot(rejections, col = colour_vector, ylab = "rejection (scaled)")
## Plot the blob wise metrics
par(mfrow = c(3, 1))
boxplot(angles %% 90, col = colour_vector, ylab = "angle (modulo 90)")
boxplot(abs(projections), col = colour_vector, ylab = "projection (scaled - absolute)")
boxplot(abs(rejections), col = colour_vector, ylab = "rejection (scaled - absolute)")
For the tip wise approach, the method is pretty straightforward from what I’ve done previously: we just measure the projection/rejection (elaboration/exploration) on each for each axis that we’ve drawn above. We thus end up with 100 projection/rejection score for each element (tip) in the tree for:
library(dispRity)
## Wrapping function
lapply.projections <- function(one_axis, trait_space) {
return(list(
projection = projections(as.matrix(trait_space), point1 = one_axis[1,], point2 = one_axis[2,]),
rejection = projections(as.matrix(trait_space), point1 = one_axis[1,], point2 = one_axis[2,], measure = "distance")))
}
## Plot rejection/projection profiles
plot.proj.rej <- function(point_list, cols, data, main) {
## Colour vector
colour_vector <- cols[morphdat$clade]
colour_vector <- adjustcolor(colour_vector, alpha.f = 1/10)
## lims
lim <- range(unlist(point_list))
plot(NULL, xlim = lim, ylim = lim, xlab = "projection", ylab = "rejection", main = main)
## Add each points
plot.points <- function(one_set, colour_vector) {
points_args <- list()
points_args$pch <- 19
points_args$col <- colour_vector
points_args$x <- one_set$projection
points_args$y <- one_set$rejection
do.call(points, points_args)
}
lapply(point_list, plot.points, colour_vector)
return(invisible())
}
## Getting the results
phylo_results_uncentered <- lapply(uncentered_major_axes$phylogeny, lapply.projections, trait_space = morphdat[, c(1,2,3)])
## Plot results
plot.proj.rej(phylo_results_uncentered, cols = colour_vector, data = morphdat, main = "Phylo elaboration/rejection")
But maybe it’s easier to boxplot that to get the info by clade
boxplot.proj.rej <- function(point_list, data, cols) {
## Get the classifier
classifier <- data$clade
sort.clade.rej <- function(data, classifier) {
return(data.frame("points" = data$rejection,
"classifier" = classifier))
}
sort.clade.proj <- function(data, classifier) {
return(data.frame("points" = data$projection,
"classifier" = classifier))
}
## Sort the data per clade
rejections <- do.call(rbind, lapply(point_list, sort.clade.rej, classifier))
projections <- do.call(rbind, lapply(point_list, sort.clade.proj, classifier))
par(mfrow = c(2,1))
boxplot(points ~ classifier, data = projections, col = cols, main = "Phylo projections\n(projected on the whole phylo effect)")
boxplot(points ~ classifier, data = rejections, rejections, col = cols, main = "Phylo rejections")
}
boxplot.proj.rej(phylo_results_uncentered, cols = colour_vector, data = morphdat)
And now we can see the differences for the within clades (i.e. the elaboration/exploration within each clade). The idea here is to distinguish between extra special taxa. For example, in mammals the platypus and the naked mole rat are special, however, the platypus is actually not special among platypuses whereas the naked mole rat is special among rodents!
## The within clades explo/elabo
gulls_results_uncentered <- lapply(uncentered_major_axes$cladegulls, lapply.projections, trait_space = morphdat[which(morphdat$clade == "gulls"), c(1,2,3)])
plover_results_uncentered <- lapply(uncentered_major_axes$cladeplovers, lapply.projections, trait_space = morphdat[which(morphdat$clade == "plovers"), c(1,2,3)])
sandpipers_results_uncentered <- lapply(uncentered_major_axes$cladesandpipers, lapply.projections, trait_space = morphdat[which(morphdat$clade == "sandpipers"), c(1,2,3)])
projection_list <- list("gulls" = unlist(lapply(gulls_results_uncentered, `[[`, "projection")),
"plovers" = unlist(lapply(plover_results_uncentered, `[[`, "projection")),
"sandpipers" = unlist(lapply(sandpipers_results_uncentered, `[[`, "projection")))
rejection_list <- list("gulls" = unlist(lapply(gulls_results_uncentered, `[[`, "rejection")),
"plovers" = unlist(lapply(plover_results_uncentered, `[[`, "rejection")),
"sandpipers" = unlist(lapply(sandpipers_results_uncentered, `[[`, "rejection")))
par(mfrow = c(2,1))
boxplot(projection_list, main = "projection within clades\n(projected on the clade specific effect)", col = colour_vector)
boxplot(rejection_list, main = "rejection within clades", col = colour_vector)
They seem a bit off as well. Maybe it’ll be worth centering their axis on their trait space centroid.
TODO.
Use the model7 with two nested clade levels?